import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import weibull_min
from scipy.optimize import curve_fit

# Datos
df1 = pd.read_csv("050525a100525.csv", sep=";", decimal=",")
df2 = pd.read_csv("171125a221125.csv", sep=";", decimal=",")

v1   = df1["velocidad"].dropna().values
dir1 = df1["direccion"].dropna().values
E1   = df1["energia"].dropna().values

v2   = df2["velocidad"].dropna().values
dir2 = df2["direccion"].dropna().values
E2   = df2["energia"].dropna().values


# Ajuste Weibull
def gamma_func(n):
    from scipy.special import gamma
    return gamma(n)

def weibull_fit_momentos(v):
    media = np.mean(v)
    sigma = np.std(v, ddof=1)
    cv    = sigma / media
    from scipy.optimize import brentq
    from scipy.special import gamma

    def cv_weibull(k):
        return np.sqrt(gamma(1 + 2/k) / gamma(1 + 1/k)**2 - 1) - cv

    k   = brentq(cv_weibull, 0.5, 20)
    lam = media / gamma(1 + 1/k)
    return k, lam

def estadisticas_weibull(k, lam, rho=1.2):
    from scipy.special import gamma
    media   = lam * gamma(1 + 1/k)
    v3_mean = lam**3 * gamma(1 + 3/k)
    v_ef    = v3_mean**(1/3)
    mediana = lam * np.log(2)**(1/k)
    moda    = lam * ((k - 1)/k)**(1/k) if k > 1 else 0
    sigma   = lam * np.sqrt(gamma(1 + 2/k) - gamma(1 + 1/k)**2)
    dpot    = 0.5 * rho * v3_mean
    return dict(media=media, v_ef=v_ef, mediana=mediana,
                moda=moda, sigma=sigma, dpot=dpot)


# Cálculos (05/05/25 - 10/05/25)
k1, lam1 = weibull_fit_momentos(v1)
stats1    = estadisticas_weibull(k1, lam1)

media1    = np.mean(v1)
sigma1    = np.std(v1, ddof=1)
TI1       = sigma1 / media1                          # Intensidad de turbulencia
v_ef1     = np.mean(v1**3)**(1/3)                   # Velocidad eficaz directa
dpot1     = 0.5 * 1.2 * np.mean(v1**3)              # Densidad de potencia directa
mediana1  = np.median(v1)

print("=== PERIODO (05/05/25 - 10/05/25) ===")
print(f"  k        = {k1:.4f}")
print(f"  lambda   = {lam1:.4f} m/s")
print(f"  Media    = {media1:.2f} m/s")
print(f"  Sigma    = {sigma1:.2f} m/s")
print(f"  IT       = {TI1:.3f}  ({TI1*100:.1f}%)")
print(f"  v_eficaz = {v_ef1:.2f} m/s")
print(f"  Mediana  = {mediana1:.2f} m/s")
print(f"  Dens.pot = {dpot1:.1f} W/m²")
print(f"  --- Weibull ---")
print(f"  Media W  = {stats1['media']:.2f} m/s")
print(f"  v_ef  W  = {stats1['v_ef']:.2f} m/s")
print(f"  Mediana W= {stats1['mediana']:.2f} m/s")
print(f"  Moda  W  = {stats1['moda']:.2f} m/s")
print(f"  Dpot  W  = {stats1['dpot']:.1f} W/m²\n")


# Cálculos (17/11/25 - 22/11/25)
k2, lam2 = weibull_fit_momentos(v2)
stats2    = estadisticas_weibull(k2, lam2)

media2    = np.mean(v2)
sigma2    = np.std(v2, ddof=1)
TI2       = sigma2 / media2
v_ef2     = np.mean(v2**3)**(1/3)
dpot2     = 0.5 * 1.2 * np.mean(v2**3)
mediana2  = np.median(v2)

print("\n=== PERIODO (17/11/25 - 22/11/25) ===")
print(f"  k        = {k2:.4f}")
print(f"  lambda   = {lam2:.4f} m/s")
print(f"  Media    = {media2:.2f} m/s")
print(f"  Sigma    = {sigma2:.2f} m/s")
print(f"  IT       = {TI2:.3f}  ({TI2*100:.1f}%)")
print(f"  v_eficaz = {v_ef2:.2f} m/s")
print(f"  Mediana  = {mediana2:.2f} m/s")
print(f"  Dens.pot = {dpot2:.1f} W/m²")
print(f"  --- Weibull ---")
print(f"  Media W  = {stats2['media']:.2f} m/s")
print(f"  v_ef  W  = {stats2['v_ef']:.2f} m/s")
print(f"  Mediana W= {stats2['mediana']:.2f} m/s")
print(f"  Moda  W  = {stats2['moda']:.2f} m/s")
print(f"  Dpot  W  = {stats2['dpot']:.1f} W/m²")


# Función histograma a intervalos de 1 m/s
def hacer_bins(v):
    vmax  = int(np.ceil(max(v))) + 1
    bins  = np.arange(0, vmax + 1, 1)
    counts, edges = np.histogram(v, bins=bins)
    centros = (edges[:-1] + edges[1:]) / 2
    freq    = counts / counts.sum()
    cum     = np.cumsum(freq)
    dur     = 1 - cum
    return centros, freq, cum, dur, edges

c1, f1, cum1, dur1, e1 = hacer_bins(v1)
c2, f2, cum2, dur2, e2 = hacer_bins(v2)


# Figura Histograma + Weibull (05/05/25 - 10/05/25)
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
fig.suptitle("05/05/25 a 10/05/25", fontsize=13)

x1 = np.linspace(0.01, max(v1) + 1, 300)
pdf1 = (k1/lam1) * (x1/lam1)**(k1-1) * np.exp(-(x1/lam1)**k1)

axes[0].bar(c1, f1, width=0.8, color="#5B9BD5", alpha=0.75, label="Datos")
axes[0].plot(x1, pdf1, color="#E07B39", lw=2, label=f"Weibull k={k1:.2f}, λ={lam1:.2f}")
axes[0].set_xlabel("v (m/s)"); axes[0].set_ylabel("f(v)")
axes[0].set_title("Histograma y ajuste de Weibull"); axes[0].legend(fontsize=8)

axes[1].bar(c1, cum1, width=0.8, color="#5B9BD5", alpha=0.75)
axes[1].set_xlabel("v (m/s)"); axes[1].set_ylabel("F(v)")
axes[1].set_title("Frecuencia relativa acumulada")

axes[2].bar(c1, dur1, width=0.8, color="#70AD47", alpha=0.75)
axes[2].set_xlabel("v (m/s)"); axes[2].set_ylabel("Fracción tiempo > v")
axes[2].set_title("Curva de duración")

plt.tight_layout(); plt.savefig("periodo1_histogramas.png", dpi=150); plt.show()


# Figura Histograma + Weibull (17/11/25 - 22/11/25)
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
fig.suptitle("17/11/25 a 22/11/25", fontsize=13)

x2 = np.linspace(0.01, max(v2) + 1, 300)
pdf2 = (k2/lam2) * (x2/lam2)**(k2-1) * np.exp(-(x2/lam2)**k2)

axes[0].bar(c2, f2, width=0.8, color="#5B9BD5", alpha=0.75, label="Datos")
axes[0].plot(x2, pdf2, color="#E07B39", lw=2, label=f"Weibull k={k2:.2f}, λ={lam2:.2f}")
axes[0].set_xlabel("v (m/s)"); axes[0].set_ylabel("f(v)")
axes[0].set_title("Histograma y ajuste de Weibull"); axes[0].legend(fontsize=8)

axes[1].bar(c2, cum2, width=0.8, color="#5B9BD5", alpha=0.75)
axes[1].set_xlabel("v (m/s)"); axes[1].set_ylabel("F(v)")
axes[1].set_title("Frecuencia relativa acumulada")

axes[2].bar(c2, dur2, width=0.8, color="#70AD47", alpha=0.75)
axes[2].set_xlabel("v (m/s)"); axes[2].set_ylabel("Fracción tiempo > v")
axes[2].set_title("Curva de duración")

plt.tight_layout(); plt.savefig("periodo2_histogramas.png", dpi=150); plt.show()


# Figura Rosa de los Vientos
def rosa_vientos(direcciones, titulo):
    bins_dir = np.arange(0, 361, 22.5)          # 16 sectores
    sectores, _ = np.histogram(direcciones, bins=bins_dir)
    sectores = np.append(sectores, sectores[0])  # cerrar el círculo
    angulos  = np.radians(np.arange(0, 360, 22.5))
    angulos  = np.append(angulos, angulos[0])

    fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(5, 5))
    ax.set_theta_zero_location("N")
    ax.set_theta_direction(-1)
    ax.bar(angulos[:-1], sectores[:-1], width=np.radians(22.5),
           color="#5B9BD5", alpha=0.75, edgecolor="white", linewidth=0.5)
    ax.set_title(titulo, pad=15)
    plt.tight_layout()
    return fig

fig3 = rosa_vientos(dir1, "Rosa de los vientos (05/05/25-10/05/25)")
fig3.savefig("rosa_periodo1.png", dpi=150); plt.show()

fig4 = rosa_vientos(dir2, "Rosa de los vientos (17/11/25 - 22/11/25)")
fig4.savefig("rosa_periodo2.png", dpi=150); plt.show()


# Figura Comparación densidad de potencia
periodos = ["(05/05/25 - 10/05/25)", "(17/11/25 - 22/11/25)"]
dpots    = [dpot1, dpot2]
colores  = ["#5B9BD5", "#E07B39"]

fig, ax = plt.subplots(figsize=(5, 4))
bars = ax.bar(periodos, dpots, color=colores, alpha=0.85, edgecolor="white")
for bar, val in zip(bars, dpots):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
            f"{val:.1f} W/m²", ha="center", va="bottom", fontsize=10)
ax.set_ylabel("Densidad de potencia (W/m²)")
ax.set_title("Comparación densidad de potencia eólica")
plt.tight_layout(); plt.savefig("comparacion_dpot.png", dpi=150); plt.show()